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Abstract 

A nonlinear time dependent fluid simulation model is developed that describes the evolution of 
magnetohydro dynamic waves in the presence of collisional and charge exchange interactions of a 
partially ionized plasma. The partially ionized plasma consists of electrons, ions and a significant 
number of neutral atoms. In our model, the electrons and ions are described by a single fluid 
compressible magnetohydrodynamic (MHD) model and are coupled self-consistently to the neutral 
gas, described by the compressible hydrodynamic equations. Both the plasma and neutral fluids are 
treated with different energy equations that describe thermal energy exchange processes between 
them. Based on our self-consistent model, we find that propagating Alfvenic and fast/slow modes 
grow and damp alternately through a nonlinear modulation process. The modulation appears to 
be robust and survives strong damping by the neutral component. 

PACS numbers: 96.25.Qr, 96.50.Dj,96.50.Bh,96.50.Tf,96.50.Xy,94.20.wf,94.30.cq 
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I. INTRODUCTION 



Alfven, fast and slow mode waves are characteristics of magnetohydrodynamic (MHD) 
plasma. While the electromagnetic fluctuations in Alfven waves propagate predominantly 
along an ambient or guide constant magnetic field, the fast/slow modes propagates isotrop- 
ically. Although the linear properties of these waves are fairly well understood 16Ml8|. 
their nonlinear evolution in a complex partially ionized environment, as occurs often in 
space/astrophysical plasmas, remains relatively uninvestigated. In partially ionized (con- 
sisting of ions, electrons and a significant neutral gas presence) space and astrophysical 
plasmas, these waves can interact with the neutral gas through processes such as charge 
exchange, collisions, ionization and recombination etc. Hence the wave properties are ex- 
pected to be modified. It is however unclear how Alfven, fast and slow waves modify the 
nonlinear processes that govern plasma fluctuations in the presence of neutral dynamics. In 
the context of cosmic ray propagation, Kulsrud & Pearce (1969) noted that the interaction 
of a neutral gas and plasma can damp Alfven waves. Neutrals interacting with plasma via 
a frictional drag process results in ambipolar diffusion (Oishi & Mac Low 2006). Ambipolar 
diffusion plays a crucial role in the dynamical evolution of the near solar atmosphere, inter- 
stellar medium, and molecular clouds and star formation. Oishi & Mac Low (2006) found 
that ambipolar diffusion can set a characteristic mass scale in molecular clouds. They found 
that structures less than the ambipolar diffusion scale are present because of the propagation 
of compressive slow mode MHD waves at smaller scales. Leake et al (2005) showed that 
the lower chromosphere contains neutral atoms, the existence of which greatly increases the 
efficiency of wave damping due to collisional friction momentum transfer. They noted that 
Alfven waves with frequencies above 0.6Hz are completely damped and frequencies below 
0.01 Hz are unaffected. They undertook a quantitative comparative study of the efficiency 
of the role of (ion-neutral) collisional friction, viscous and thermal conductivity mechanisms 
in damping MHD waves in different parts of the solar atmosphere. They found that a cor- 
rect description of MHD wave damping requires the consideration of all energy dissipation 
mechanisms through the inclusion of the appropriate terms in the generalized Ohms law, the 
momentum, energy and induction equations. Padoan et al (2000) calculated frictional heat- 
ing by ion-neutral (or ambipolar) drift in turbulent magnetized molecular clouds and showed 
that the ambipolar heating rate per unit volume depends on field strength for constant rms 
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Mach number of the flow, and on the Alfvenic Mach number. 

The role of ion-neutral collisions has been investigated extensively by Balsara (1996) in 
the context of molecular clouds. Here the momentum equation for the plasma is governed 
predominantly by the slow comoving massive neutrals that tend to dissipate the Cloud's 
magnetic field. Because of the larger neutral to ion mass ratio, ion momentum is entirely 
dominated by the neutral drag. Hence the inertial terms were ignored in plasma momentum 
equation. Balsara (1996) found that slow waves propagate without significant damping on 
short wavelengths, while the fast and Alfven waves undergo rapid damping in super- Alfvenic 
regimes. 



In the heliospheric boundary regions, particularly the outer heliosheath [13|, |20j, the 
plasma is mediated by interstellar neutral gas by virtue of charge exchange. For example 
in the local interstellar medium (LISM), the low density plasma and neutral hydrogen (H) 
gas are coupled primarily through the process of charge exchange. In this process, both the 
total (proton+atom) momentum and energy are conserved. On sufficiently large temporal 
and spatial scales, a partially ionized plasma is typically regarded as equilibrated; this is the 
case for the LISM, but not for the outer heliosheath, for example. 

It is worth mentioning that the charge exchange process in its simpler (and leading order) 
form can be treated like a friction or viscous drag term in the fluid momentum equation, 
describing the relative difference in the ion and neutral fluid velocities. The drag imparted 
in this manner by a collision between ion and neutral also causes ambipolar diffusion, a 
mechanism used to describe Alfven wave damping by cosmic rays (Kulsrud & Pearce 1969) 
and also discussed by Oishi & Mac Low (2006) in the context of molecular clouds. 

In this Letter, we focus on understanding the propagation characteristics of MHD waves 
by considering collisional and charge exchange interactions thus including both simultane- 
ously. Our model is more general and fully nonlinear, compared to that of Balsara (1996) 
which ignores nonlinear intertial terms in the plasma momentum equation. We retain this 
term and carry out fully time-dependent simulations in two dimensions. We treat plasma 
and neutral fluid in a self-consistent manner by employing the complete nonlinear time- 
dependent fluid equations for both the fluids. In Section 2, we discuss the equations of a 
coupled plasma-neutral model, their validity, the underlying assumptions and the normaliza- 
tions. Section 3 describes results of our nonlinear, coupled, self-consistent fluid simulations. 
We find that the propagation characteristics of Alfven waves are altered signifincantly by 
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the combined action of charge exchange and collisional interaction processes. The latter 
leads to a modulated growth and damping of the Alfven and fast/slow modes. In Section 
4, we summarize our results and describe physical processes that may be responsible for the 
modulation process. 

II. PARTIALLY IONIZED MODEL EQUATIONS 

We assume that fluctuations in the plasma and neutral fluids are isotropic, homogeneous, 
thermally equilibrated and turbulent. A constant mean magnetic field is present. Local 
mean flows and nonlinear structures may subsequently be generated by self-consistently ex- 
cited nonlinear instabilities. The boundary conditions are periodic. This allows us to assume 
an infinite partially ionized plasma where the length scale of characteristic fluctuations (£) 
is much smaller than the size of plasma or computational box (L) i.e. I 1 < L. We further 
assume that plasma particles interact with the neutral gas via collisions as well as charge 
exchange. The intrinsic magnetized waves supported by the plasma scatter plasma parti- 
cles in a random manner thus sufficiently isotropizing the plasma and neutral distribution 
functions [8J . The latter enables us to use fluid descriptions for both the plasma and neutral 
gas. Most of these assumptions are appropriate to realistic space and astrophysical turbu- 
lent flows. They allow us to use MHD and hydrodynamic descriptions for the plasma and 
the neutral components respectively. In the context of the LISM and the outer heliosheath 
plasmas, the plasma and neutral fluid remain close to thermal equilibirium and behave as 
Maxwellian fluids. Our model simulates a plasma-neutral fluid that is coupled via collisions 
and charge exchange in astrophysical plasmas. The fluid model describing nonlinear turbu- 
lent processes, in the presence of charge exchange and collision, can be described in terms 
of the plasma density velocity (U p ), magnetic field (B), pressure (P p ) components 

according to the conservative form 

dF 

+ V • Q p = Q p , n , (1) 
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Note the presence of the source terms Q that couple the plasma self-consistently to the 
neutral gas. The above set of plasma equations is supplimented by V ■ B = and is coupled 
self-consistently to the neutral density (p n ), velocity (V n ) and pressure (P n ) through a set 
of hydrodynamic fluid equations, 
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Equations (JT[) to ([2]) form an entirely self-consistent description of the coupled plasma-neutral 
turbulent fluid in a partially ionized medium. 

Several points are worth noting. The charge-exchange momentum sources in the 
plasma and the neutral fluids, i.e. Eqs. (JTJ) & (j2J), are described respectively by 
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terms QM, P , n (U P|1 Vn, p P , Pn, T n , T p ) and QM,n, P (V n , U p , p p , p n , T n , T p ). These expressions, 



described in 
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131 ] . have the following form. 
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A swapping of the plasma and the neutral fluid velocities in this representation corresponds, 
for instance, to momentum changes (i.e. gain or loss) in the plasma fluid as a result of charge 
exchange with the neutral atoms (i.e. QM,p,n(U p , V n , p p , p n , T n , T p ) in Eq. (JT])). Similarly, 
momentum change in the neutral fluid by virtue of charge exchange with the plasma ions is 
described by Q M ,n, P (Vn, U p , p p , p n , T n , T p ) in Eq. (jSJ). 

The plasma-neutral collisional forces are modeled by employing the following forcing term 
in Eqs. © & ©. 

Fp, n = z/ppp n (U p - V n ) = -F n>p , (5) 

where v is the collision frequency. The collisions modify plasma and neutral momentum and 
energy, but not the density. The negative sign before F n p in Eq. (jSJ) ensures conservation 
of total momentum. We derive modified equations of energy for both fluids following the 



treatment described in 



12] . The collisional force is self-consistently calculated in our model. 



Note also the charge exchange source term in momentum and energy of both the fluids. In 
the absence of collision and charge exchange interactions, the plasma and the neutral fluid are 
de-coupled trivially and behave as ideal fluids. Like collisions, charge-exchange interactions 
modify the momentum and the energy of plasma and the neutral fluids and conserve density 
in both the fluids (since we neglect photoionization and recombination). Nonetheless, the 
volume integrated energy and the density of the entire coupled system will remain conserved 
in a statistical manner. The conservation processes can however be altered in the presence 
of any external forces. These can include large-scale random driving of turbulence due to 
external forces such as supernova explosions, stellar winds, etc and instabilities. Finally, 
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the magnetic field evolution is governed by the usual induction equation, i.e. Eq. ((T]), that 
obeys the frozen-in-field theorem unless some nonlinear dissipative mechanism introduces 
small-scale damping. 

The underlying coupled fluid model can be non-dimensionalized straightforwardly using 
a typical scale- length (£ ), density (po) an d velocity (vq). The normalized plasma density, 
velocity, energy and the magnetic field are respectively; p p = p p /po, U p = XJ p /vq,P p = 
Pp/pofg,B = B I 'vq-s/pq. The corresponding neutral fluid quantities are p n = p n /p ,XJ n = 
XJ n /vQ,P n = P u /pqVq. The momentum and the energy charge-exchange terms, in the nor- 
malized form, are respectively Q m = Qm^o/po^O' Q e = Qe^o/Po v o- The non-dimensional 
temporal and spatial length-scales are t = tvo/£o,x = x/%- Note that we have removed 
bars from the set of normalized coupled model equations (CQ) & (T5]). The charge-exchange 
cross-section parameter (a), which does not appear directly in the above set of equations, 
is normalized as a = n £ a, where the factor Uq£q has dimension of (area) -1 . By defining 
rio,£ Q through a ce = l/n £o = k^ e , we see that there exists a charge exchange mode (k ce ) 
associated with the coupled plasma-neutral turbulent system. For a characteristic density, 
this corresponds physically to an area defined by the charge exchange mode being equal to 
(mpf) 2 . Thus the larger the area, the higher is the probability of charge exchange between 
plasma ions and neutral atoms. Therefore, the probability that charge exchange can directly 
modify those modes satisfying k < k ce is high compared to modes satisfying k > k ce . Since 
the charge exchange length-scales are much smaller than the turbulent correlation scales, 
this further allows many turbulent interactions amongst the nonlinear turbulent modes be- 
fore they undergo at least one charge exchange. An exact quantitive form of sources due to 
charge exchange in our model is taken from Shaikh & Zank (2008). 



III. EVOLUTION OF MHD WAVES 



We now investigate the nonlinear evolution of the coupled plasma-neutral system de- 
scribed by Eqs. (JTJ) & ([2]). Our goal is to understand the evolution of MHD waves in a 
partially ionized turbulent plasma, mediated by the complex coupling between the two dis- 
tinguishable fluids. For this purpose, we initialize the field variables (i.e. velocity, magnetic, 
density, energy and pressure fields) in the plasma and neutral fluids with a random, un- 
correlated, out-of-phase distribution in space. The initial ratio of the plasma and neutral 



7 



density is unity. The plasma and neutral fluids initialized in this manner characterize spa- 
tially local turbulence in the partially ionized astrophysical plasma. These random fields 
develop spatially and temporally according to the set of Eqs. (JTJ) &; (ED and forces due to 
charge exchange and collision dictate nonlinear interaction processes. In ideal MHD and 
hydrodynamic fluids (i.e. without the coupling forces in Eqs. (JTJ) & known quadratic 
conserved quantities like energy, helicity, vorticity etc predominantly govern the spectral 



^ange. Their role is 
5|, ISJ. By contrast, 



transfer of energy across disparate characteristic scales in the inertia 
well studied within the paradigm of statistical theories of turbulence 
the complex nature of the coupling forces in Eqs. (jTJ) & (j2J) poses a formidable hurdle to 
deriving quadratic conserved quantities. It is not clear to us as how these quantities will 
influence spectral transfer in the inertial range. Furthermore, the driving forces in Eqs. ([T]) 
& (j2J) may presumably be operative at any length and time scales. It is unclear how this 
will alter the conventional dissipative and diffusive processes in a partially ionized plasma. 
Unlike the constant drag force (i.e. not changing in space and time) used before for a neu- 
tral fluid; and which acts to damp Alfven waves in the partially ionized plasma (Kulsrud & 
Pearce 1969, Balsara 1996, Shaikh 2010), it is not clear how the complex, nonlinear, time 
and space dependent driving forces in Eqs. (JTJ & (j2J) will modify the propagation and inter- 
actions characteristic of Alfven, fast and slow MHD waves or even compressive sound waves 
in the coupled neutral gas. 

To address the above issues, we have developed a two-dimensional (2D) nonlinear fluid 
code to numerically integrate Eqs. (0Q) & (j2J). The 2D simulations are not only compu- 
tationally less expensive (compared to a fully 3D calculation), but they offer significantly 
higher resolution (to compute inertial range turbulence spectra) even on moderately-sized 
clusters like our local Beowulf system. The spatial discretization in our code uses a discrete 
Fourier representation of turbulent fluctuations based on a spectral method, while we use 
a 4th-order Runge-Kutta method for the temporal integration. All the fluctuations are ini- 
tialized isotropically with random phases and amplitudes in Fourier space. A mean constant 
magnetic field Bq is assumed along the y-direction. Our algorithm ensures conservation of 
total energy and mean fluid density per unit time in the absence of charge exchange and 
external random forcing. Additionally, V • B = is satisfied at each time step. Our code 
is massively parallelized using Message Passing Interface (MPI) libraries to facilitate higher 
resolution. The initial isotropic turbulent spectrum of fluctuations is chosen to be close to 
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FIG. 1: Evolution of modal energy in Alfven and fast/slow waves is shown by their respective 
average mode behavior. The oscillatory behavior of the average mode associated with Alfven and 
fast / slow indicates that charge exchange and collisions with neutrals damp the Alfven and fast / slow 
wave which are then reexcited by an instability (growth). This process repeats itself during the 
entire evolution. 




time 



FIG. 2: Energy exchange rate between the plasma and neutral fluids shows that the rate of transfer 
of energy between the two fluids is constant. We plot the corresponding energy equations from 
Eqs (1) & (2). 

k~ 2 with random phases in all three directions. The choice of such a (or even a flatter than 
-2) spectrum does not influence the dynamical evolution as the final state in our simulations 
progresses towards fully developed turbulence. 

From our time-dependent nonlinear fluid simulations, we determine the quantative evolu- 
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tion of Alfven and fast/slow waves in a partially ionized environment. For this purpose, we 
distinguish Alfvenic and non- Alfvenic, i.e. corresponding to the compressional or due to slow 
and fast magnetosonic modes, contributions to the turbulent velocity fluctuations. Thus, 
we introduce diagnostics that distinguish between the modes corresponding to Alfvenic and 
slow/fast magnetosonic fluctuations. Since the Alfvenic fluctuations are transverse 
the propagation wave vector is orthogonal to the velocity field fluctuations i.e. klU, and 
the average s pec tral energy contained in these (shear Alfvenic modes) fluctuations can be 



computed as [14 1 



The above relationship yields a finite spectral contribution from the |k x Uk| characteristic 
turbulent Alfvenic modes. It is essentially a modal energy associated with the Alfvenic 
fluctuations that obey Eq. (J5]) in the inertial range turbulent spectrum. On the other hand, 
fast/slow (i.e. compressive) magnetosonic modes propagate longitudinally along the velocity 
field fluctuations, i.e. k || U and thus carry a finite component of energy corresponding 
only to the zk • Uk part of the velocity field, which can be determined from the following 
relationship 

The expression for k c essentially describes the modal energy contained in the non-solenoidal 
component of the MHD turbulent modes. 

The evolution of an average mode associated with Alfvenic and compressive waves is 
shown in Fig. (JTJ). We vary charge exchange strength and collision parameter in our sim- 
ulation to examine their effects on the propagation of Alfvenic and compressive modes in 
a partially ionized plasma. Clearly, the two processes operate on different time and length 
scales and are self-consistently modeled in our simulations. We find that charge exchange 
and collisional interactions jointly modulate Alfven and fast/slow modes. This is shown in 
Fig. (JTj) for 512 2 modes in a two dimensional box of length 2ir x 2n. The other parameters 
in our simulations are charge exchange wavenumber k/k ce ~ 0.01, fixed time step dt = 10~ 3 , 
and collision parameter v ~ 0.001. The background normalized constant magnetic field 
Bq = 0.5. Our simulations are fully nonlinear and the ratio of the mean and fluctuating 
magnetic fields 5B/B ~ 1. 

The modulation process depicted in Fig. (JT]) can be understood on the basis of charge 

10 



exchange and collisional interactions, together with convective nonlinearities. During the 
modulation, there exists a rapid onset of the average mode associated with Alfvenic waves. 
This onset is triggered essentially by linear instability process. Our linear stability analysis 
(not described here) indicates that the underlying coupled plasma-neutral system possesses 
several unstable modes. These modes account for the growing as well as damping of Alfvenic 
and fast/slow compressive waves. As illustrated in Fig. (JTJ, a concurrent damping of the 
fast/slow compressive mode preceding the growth of Alfvenic mode occurs. When the two 
modes reach their extremal rise or fall, they reverse their behavior. The fast/slow compres- 
sive mode begins to rise at the expense of damping of Alfvenic mode. This process continues 
non periodically and repeats itself in time thereby exhibiting a nearly oscillatory behavior. 
It follows from Fig. ([rj that the two modes, i.e. Alfvenic and fast/slow compressive, reg- 
ulate each other in a predator-prey manner. Dynamically, the Alfvenic mode grows at the 
expense of fast/slow mode. When Alfvenic mode reaches its maximum amplitude, where 
the fast/slow mode remains at its lowest magnitude, the latter eats up the Alfvenic mode 
and vice versa. This process is repetitive. It is also noteworthy that the overall amplitude 
of both the modes decays in time and regrows later. 

The growth and damping of wave energy can further be understood from the energy 
transfer rates (dE/dt) between the plasma and neutral fluids. This is shown in Fig. ([2]). It 
is evident from the figure that the total energy of the plasma increases at the expense of the 
neutral. When the waves damp, the total energy of the plasma decreases while it increases 
for the neutrals. The energy transfer rates also exhibit a periodic behavior following the 
modulation of energy in Alfven and fast / slow modes. Nonetheless, the energy transfer rates 
try to approach nearly zero i.e. dE/dt — > 0, but are not entirely zero due to the modulation 
process. This is an interesting result that differs from pure (i.e. uncoupled) plasma or 
neutral fluid turbulence {2, 8]. The latter exhibits dE/dt ~ thereby leading to a constant 
transfer of energy across inertial range fluctuations. This further demonstrates that the total 
energy of the neutral and plasma fluid is nearly constant during the steady state (See Fig 2). 
What we find in our simulations is a modulated-energy-transfer process that is completely 
different from the constant energy transfer rate. The energy transfer rates in Fig. (T5]) thus 
provide a self-consistent description of the modulation process observed in our simulations. 
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IV. CONCLUSION 



In summary, we developed a two dimensional self-consistent model of plasma and neu- 
tral fluids that is coupled through charge exchange and collisional interactions. Notably, 
we find that these interactions not only modify the convective nonlinear interactions in the 
coupled Eqs. (fll) &: ([2"]), but they also influence the propagation characteristics of MHD 
waves significantly. To determine the evolution of characteristic MHD waves, we distinguish 
the Alfvenic and non-Alfvenic contributions of MHD modes in the partially ionized turbu- 
lent environment. It is found that charge exchange and collisional interactions lead to a 
modulation of Alfvenic and fast/slow compressive modes. This is unlike purely collisional 
processes described by fl and A that tend to damp Alfvenic modes. By contrast, inclusion 
of the charge exchange interactions leads to an entirely different scenario and it alters the 
propagation characteristics of MHD waves. Charge exchange, in the presence of collisions 
and convective nonlinear interactions, leads to an alternate growth and damping (i.e. modu- 
lation) of Alfvenic and fast / slow compressive modes. It is further noted that the modulation 
of waves is mediated by the conservation of momentum because of the characteristic wave 
velocity that is involved in Eqs. (0) & O- It is the momentum transfer function (Q m ) that 
is responsible for regulating the rise/fall of the amplitude of the two waves and maintains 
the conservation process. This point is further consistent with Fig (1). We finally comment 
on the conservation of momentum and energy in our simulations. It seems from Fig (2) 
that the term dE/dt hovers around zero, but does not fully converge to zero. We learn that 
due to the small scale self-consistent dissipation (or numerical dissipation), it is difficult 
to observe a 100% conservation. We continue to work to improve this conservation in our 
model. 

Our results should find application to a variety of astrophysical environments in which a 
partially ionized plasma is typical. Examples include the outer heliosheath formed by the 
interactions of the solar wind with the local interstellar medium, the magnetic collapse of 
molecular clouds and star formation and t he g eneral transfer of energy in partially ionized 



plasmas surrounding other stellar systems 



19]. 
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